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Abstract 

We investigate the macroscopic quantum tunneling of the attractive Bose-Einstein con- 
densate. Within the effective Lagrangian framework, we find bounce solutions and explicitly 
^ | calculate the decay rate of the condensate trapped in a cylindrically symmetric potential. In 

particular, in the case where the number of condensed bosons is slightly below a certain critical 
number, we present a detailed analysis of the bounce solutions and discuss the approximations 
employed in our calculations. The effects of finite temperatures and the shape of the trapping 
potential are evaluated. 

e. 

ii 1 Introduction 

c 

q \ Macroscopic quantum tunneling is an interesting subject in many areas of physical sciences including 
low-temperature physics, atomic physics and nuclear physics. Recent realization of the Bose-Einstein 
condensate of trapped alkali atoms may provide a good testing ground for the investigation of this 
problem [l]. 

In this paper we will discuss the macroscopic quantum tunneling of the Bose-Einstein condensate 
with attractive interactions. The dynamics of the condensate is successfully described by the Gross- 
Pitaevskii (GP) equation 0, [|. The s-wave scattering length a entering the GP equation can be 
positive or negative, its sign and magnitude depending crucially on the details of the atom-atom 
interaction. In the case of 7 Li, the interaction is attractive and the scattering length is known to be 
a = —1.45 ± 0.04 nm |4], ||. The attractive interaction causes the condensate to collapse upon itself. 
When the trapping potential is included, however, the destabilizing influence of the interaction is 
balanced by the zero-point kinetic energy, thereby allowing a metastable condensate to form || |7], §]. 
Perez-Garcia et al. |5| have investigated the GP equation by using a time-dependent variational 
ansatz for the condensate wave function. Their results reproduce quite accurately the low energy 
excitation spectrum of the condensate obtained by numerical simulations of the GP equation. We 
will apply this variational technique to the macroscopic tunneling of the metastable condensate of 
7 Li. When the trapping is spherically symmetric, Ueda and Leggette have evaluated the tunneling 
decay rate at zero temperature (see also [J7|, [TO]). In this paper we develop their analysis and 
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explicitly write down the decay rate in the case of a cylindrically symmetric trapping potential and 
further finite temperatures. 

In Section 2, according to ||, we derive an effective Lagrangian describing the Bose-Einstein 
condensate of 7 Li, and summarize the data of the ground state energy that we shall need in the 
calculations of tunneling. In Section 3, we present a detailed analysis of bounce solutions. Using the 
effective Lagrangian and with the help of numerical simulations, we find the bounce solutions. We 
next consider the special situation, where the number of condensed bosons is slightly below a certain 
critical number. Then the effective Lagrangian reduces to a simple one-dimensional Lagrangian by 
appropriate approximations. We present an analytic solution for the bounce within this situation, 
and explicitly calculate the decay rate of the metastable condensate. We also evaluate the decay rate 
at finite temperatures and predict a critical temperature, where the rate crosses over from quantum 
tunneling to thermal hopping. Section 4 is devoted to the summary of our findings. 



2 Model 

We consider gases of 7 Li atoms trapped in a cylindrically symmetric harmonic potential 

V(x,y,z) = ^mu 2 {x 2 + y 2 + X 2 z 2 ), (1) 

where A represents the asymmetry parameter of the trapping potential. The dynamics of the con- 
densate is described by the GP Lagrangian 

ih ( dip* dip , A h 2 2 2 2nh 2 a. , l4 

In order to obtain the evolution of the condensate wave function, we assume the Gaussian form 
for the wave function according to 0: 



ip(x,y,z,t) = A(t) J] exp 



a=x,y,z 



2W a {t) 2 
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ix a a a (t) + ix a /3 a (t) 



(3) 



This trial function includes the time- dependent variational parameters, ry = (r] x ,r]y,r] z ) (center co- 
ordinate), W = (W x , Wy, W z ) (width) and the phase parameters ct = (a x ,a y ,a z ), (3 = (fl x , f3 y , (3 Z ) 
which correspond to the canonically conjugate "momentums" to rj and W. The wave function ip is 

normalized by the number of condensed bosons N = J \ip\ 2 d 3 x, so that the parameter A (amplitude) 

is given by 



Substituting into (0) and further integrating the GP Lagrangian over space coordinates, one 
obtains an effective quantum mechanical Lagrangian 

^e// = J2 (PaVa + K a W a ) - TUffiVa, W a ,p a , K a ), (5) 
a=x,y,z 
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where p a and K a are the canonically conjugate momentums to rj a and W a defined by 

p a = hN(a a + 2r] a p a ), (6) 
K a = hN(3 a W a . (7) 

The Hamiltonian TC e ff = 7~Lo + H consists of two parts: the first part simply describes the 
harmonic oscillation of the center of the condensate 

^o= E ^A + ^tt + ri + W), (8) 

a=x,y,z ^" iiv ^ 

and the remaining part TC describes the evolution of the widths of the condensate 

n = E ±K* a + U(W) (9) 

a=x,y,z 

with 

r -\. mNis 2 . Tir2 - x2 „. h 2 N ( 1 1 1 \ ah 2 N 2 1 

^ W = 7— ( W x + W» + X W z + -r~ 777? + 777? + 7777 + 7= f77777T7T - 10 

4 y 4m \W£ W 2 W^ 2 / v^fmW^ 

It is convenient to introduce the scales characterizing the trapping potential : (a) length scale ao = 



U/mv, (b) energy scale e = Uv/2, (c) time scale v 1 . By using these units we define dimensionless 
quantities, £ = clq 1 ^, X = 1 W and r = vt. Then the Lagrangian (^) is rescaled as follows : 



e 1 C eff = L + L, (11) 



where 



and 



^o = iV(^j -iV(^ + ^ + A 2 e 2 2 ), (12) 



= i(* + Y> + X*) + l(± + ± + ±) + ^ (14) 

with P = ^2/irNa/ao < 0. We now focus our attention on the ground state energy of the condensed 
Bose system. Under the present analysis, the ground state energy can be calculated by finding the 
critical points of U and the eigenvalues of the Hessian matrix H a b = d 2 U / dX a dX^ evaluated on the 
critical points. 
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Critical points 



The critical points are given by the solutions to 



X = Y, 
1 P 

+ x 2 z 3 : 
1 p 



A 2 , 
= I. 



(15) 
(16) 

(17) 



The solutions are classified by the critical value P* of the parameter P; when \P\ > \P*\, there are 
no critical points. When \P\ < \P*\, there are two critical points, one stable (Morse index = 0) and 
the other unstable (Morse index = 1) 0. The critical value P* satisfies in addition to the equations 
(H) © (0) also 



1 P 2 



X 2 Z 3 2X 6 Z 4 



4X 2 



It 



which can be derived from the condition ex = (see (0)). Thus, P* and the corresponding coordinate 
X* = (X*, Y*, Z*) are uniquely determined as a function of the asymmetry parameter A. Indeed we 
have P* = — 4/5 5 / 4 , X* = 5 _1 ^ 4 (1, 1, 1) for A = 1, and general solutions are provided in Fig.|l]. It 
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Figure 1: A-dependence of P* and X* = (X*,Y*, Z*" 



should be noticed that for P — > P* the stable critical point X s and the unstable critical point X u 
take the following asymptotic forms: 



X f 



X* ± k(l - P/P*fl 2 E + Oil - P/P* 



where E = (— P 3 * 2 , — P 3 * 2 , 1) and the coefficient k is given by 



P* 1 (2P 



li 



2A 2 (1 - P 



p* 
r 2Z , 



(19) 



(20) 
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with P*. = P*/A{X*)\Z 



Eigenvalues of Hessian matrix 

For the eigenvalue problem of the Hessian matrix evaluated on the critical points 

He A = e\e A (A = T, N, B), 

we have the following results 0: 
fa) T-direction 



2 A + 1 - P 23 - j8(P 32 y + (1 - A 2 + P 23 ) 2 



e T 

with T 



A T 

i (->> 



(T,T,P; 



32/ 



A + 1 + P 23 - J8(P32) 2 + (1 - A 2 + P 23 ) 2 



(b) A^-direction 





2 (A 2 






e N = 




with N = 


K- 







2 (A 2 + 1 - P 23 + v/8(P 32 ) 2 + (1 - A 2 + P 23 ) 2 ) 
(N,N,P 32 ), 



fc) P-direction 



e 2 , = 4(1 - 2P 



41; 



e B = -^(1,-1,0). 



Here we used the notation 



P 



(\P\ < \P'\), 



and 



./'. V - -< P'->2 ) ' + 1 



(21) 



A 2 o / d ^2 



(P 23 + 1 - A ) =f (P 23 + 1 - A ) a/8(P 32 ) 2 + (1 - A 2 + P 23 ) 2 



(22) 
(23) 
(24) 

(25) 
(26) 
(27) 

(28) 
(29) 

(30) 
(31) 



by the normalization e^-e^ = 5ab- It should be noticed that the eigenvalue e T is positive (negative) 
for the stable (unstable) critical point and the other eigenvalues are all positive. These results imply 
the following ground state energy 



e l E g = NU(X S ) + (2 + A) + £ e A (X s ). 

A=T.N,B 



(32) 



The first term is the potential energy evaluated on the stable critical point X s , and the second-third 
terms represent the zero-point energy coming from collective excitations of the condensate. 
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3 Macroscopic quantum tunneling 



In this section we argue the macroscopic quantum tunneling of the Bose condensate using the La- 
grangian (|T3|). The stable critical point X s of the potential U(X) represents a metastable condensate 
since the parameter P in U(X) is negative, and so the ground state energy will have an (exponen- 
tially small) imaginary part in addition to ( |32|) if we take account of the tunneling. The decay rate 
of the metastable condensate is determined from 

T = ^Im£ g . (33) 

We will calculate the decay rate by using the WKB approximation. Since the Lagrangian (|i"3|) 
includes a macroscopic quantity iV representing the number of condensed bosons, we must be careful 
for the choice of a small parameter h controlling the validity of the WKB approximation. The precise 
value of h is given by fl52|), and the decay rate is of the form 

r~Aexp(-^), (34) 

where S c \ is the Euclidean action evaluated at bounce solution and A the square root of the deter- 
minant of the second variation around the bounce solution, with the zero - mode removed. 



Zero temperature 

We start with the Euclidean action: 



h 2 7-oo \2\dr J 



where U(X) is the potential given by ([14]) . The bounce solution is the classical solution to the 
equations of motion 

d 2 1 1 

X-X + — + P-— = 0, (36) 



dr 2 X 3 X 2 YZ 

f Y ~ Y + T* +P xk = °> < 37 > 

<]2 Z-X 2 Z+^ + P-l- = 0, (38) 



dr 2 Z 3 XYZ 2 

subject to the boundary condition 

lim X(t) = X s (stable critical point). (39) 

T — >±00 

In Fig.|2|. we show the behavior of bounce solutions obtained using numerical simulations. 

Let us investigate analytically the system (|35D by choosing a parameter P near the critical value 
P*. Then the bounce solution X(,(t) is restricted in the neighborhood of the stable critical point 
X s . Indeed, the equation (|i~9| ) gives the estimation, |Xb(r) — X s \ ~ \X S — X U \ ~ 0((1 — P/P*) 1 ^ 2 ). 
In the following text we will assume 

5 = 1 - P/P* ~ 1(T 3 , \P\<\P*\. (40) 
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(a) A = 0.5, P = -0.5 



(b) A = 2.0, P= -0.4 



Figure 2: Behavior of the bounce solution. The bold-faced curve connecting the two points, X s (stable 
critical point) and X t (turning point) corresponds to the bounce solution. The solid curves represent the 
contours of the potential U(X). Parameter values : 5 = 1 — P/P* = 0.144 for (a) and 5 = 0.135 for (b). 



This parameter region is particularly interesting; as seen later on the value is of the order 

one, in this region, though the prefactor N in the action is very large (the number of atoms used 
in the experiment at Rice University is of order 10 3 [ IT , |I~2|j). Thus we can expect to observe the 
macroscopic quantum tunneling by experiments. 

We now introduce a new coordinate x = (xt,xn,xb) around X s 

X = X S + J2 x Ae A (41) 

A=T,N,B 

and expand the potential 

U{X)=U{X s ) + \ J2 44 + E c nml x%x%x l B + ---. (42) 

1 A=T,N,B n+m+l=3 

It should be noticed that the eigenvalue approaches to zero for 8 — > 0. Indeed, we can evaluate 
the behavior of near P* using the exact formula (|22|) : 

e T = a5 1/i + 0(5 3/4 ), (43) 

where 



A 2 + 1 - P 2 3 

On the other hand, eigenvalues and e# can be approximated by ( p5|) and (p8|) evaluated on P*, 
and these values become extremely large compared with when the parameter 5 approaches to 
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zero. This means that the direction of the initial (infinitesimal) velocity of the bounce solution 
is given by the eigenfunction ey. Thus the trajectory of the bounce solution is mainly described 
by x(t) = xt{t), i.e. T-component of the coordinate x(t), and remaining components xat(t) and 
xb(t) give higher order corrections. More precisely, using fl42| ) and (|43|) , we can evaluate the bounce 
solution as xt(t) ~ 0{5 1 ^ 2 ), xn{t) ~ 0{5) and xb(t) = by the symmetry of the equations of 
motion (if we specialize to the spherically symmetric trapping potential, the iV-component xjy(r) 
exactly vanishes). We now approximate ([35]) by one-dimensional quantum mechanical action: 



S E N (\ ( dx\ 2 1 2 2 c ; i 

T ~~2 i-oo ( 2 ( <i( ) + ~i T + 3! ' ' (45) 



From (p3|) and (|2~4[) , the coefficient c(< 0) is given by 



1_ f /6» 3 C7 _ d 3 U \ „,d 3 U „„„„ d 3 U „ d 3 U 

which takes the following asymptotic form 



c " M \ 2T " Idx*- + W + Ps W + 18T2p32 dxdYdz + 6TP32 axdzi 



, (46) 

X=X S 



c = -12(P 2 * 4 + 4P 4 * 1 P 2 * 4 - 2A 2 P 4 * 2 )(1 + 2(P 3 * 2 ) 2 )- 3/2 + 0(S 1/2 ). (47) 

It is convenient to introduce new scales characterizing the quantum tunneling : according to Fig. 3 
we define 

(a) length scale R = = !^V/ 2 (1 + Oti 1 ' 2 )), (48) 

\c\ \c\ 

(b) energy scale U = ^4 = ^V 2 (l + O^' 2 )). (49) 
Then we have a natural time scale 

To = ^^ = ^ r V4 (1 + 0(5 i/ 2)); wo = ./^ i (50) 

(2U /Nm) 1/2 vol 1 1 V 2 1 ; 

representing the "tunneling time" . 
Now the action (|45T) is of the form 



S " ' r ds[lrf] + #(?)), U(q) = hjtf(l-q), (51) 



2 

7/ // 7_oo"" I 2 \ r/.s j - ""/' 2" 

where we have used rescaled quantities, q = (a /R )x and s = (1/VT )t. The prefactor (effective 
Plank constant) 

A -4^ r5/1(1+0( ' 1/2)) (52) 

is a dimensionless parameter controlling the validity of the WKB approximation . The equations of 
motion can be easily integrated yielding the well known bounce solution 

q b (s) = sech 2 (^) . (53) 




Figure 3: Potential profile for 5 — > 0. The potential U(W) given by (|io|) is approximated by a one- 
dimensional potential U(R) = eoN((e^/2)x 2 + (c/3!)x 3 ) with R = clqx. The potential U(R) has a metastable 
minimum at R = and a barrier of height Uq = U(R m ),Rm = 2aoe^/\c\. 



Using the WKB approximation [T3|], we obtain the decay rate 

To = I 4\/^ exp ( ) ) (1 + Oih))^- 1 

and 




with the bounce action S c ] = — ujq. It follows from 

15 

contribution to T Q is given by 



that for 5 



cxp 

Here, the coefficients A and -B are functions of the asymmetry parameter A: 

~Ya 7/2 „ 12a 5 




5c 2 



(54) 

the leading 
(55) 

(56) 



which can be calculated by ( |44| ) and (p]). Fig|| shows the A-dependence of these coefficients. The 
spherically symmetric trapping potential (A = 1) minimizes the function B and its value is 4.58 in 
excellent agreement with the result of ||. The functions A,B remain relatively constant for A < 1 
but they grow for A > 1. For 5 —>■ the tunneling exponent and the prefactor vanish according to 
<5 5//4 and <5 7//8 , respectively P, [TD|. We find that this scaling law is universal, independently of the 
shape of the harmonic trapping potential. 



Finite temperature 



In the case of finite temperature the bounce solution is given by a periodic solution, i.e. the 
classical solution in the potential —U(q) with energy — E (0 < E < 1). From Fig.|| the solution 
takes the form 

(57) 



qb(s) = <? 2 - (<?2 - <?i)sn 2 ( 



^0 , 



q s;m 
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with the elliptic modulus m 
of the first kind : 



<?2 - go 



and the period h(5 is given by the complete elliptic integral 



hp = A K(m), K(m) = C 

^o\/<?2 - go Jo 



dx 



1 — x 2 )(l — mx 2 ) 



(5f 



This solution reduces, of course, to the previous solution (|53|) for _E = 0. The corresponding 
bounce action is evaluated as 



Sr. 



W + hpE, 



(59) 



where 



W 
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VQ2 - <?o[2(<? + ?i + & - 9o9i - Q0Q2 - qiQ2)E{m) + (<?i - <?o)(2g - gi - g 2 )^(m)] (60) 



E(m) is the complete elliptic integral of the second kind ). The fluctuation modes about the bounce 



solution include a zero mode = (fc>(s). Then the determinant factor A in (34) is calculated from 
the Gelfand-Yaglom formula |15L |16| : 



sinh (uqs) 



ds' 



s=f3h/2 

Thus we obtain the finite temperature decay rate due to quantum tunneling : 

-S c i 



(61) 



T((3) = (A((3)exp 



h 
10 



[l + O(h))T -\ 



(62) 




where 



with 



I (gg ~ 9o) 3/4 (92 - 9i)(l - m 2 ) 
2nh (a(m)E(m) + b(m)K(m)y/ 2 



sinh 



(63) 



a(m) = 2(m 4 -m 2 + l), (64) 
b(m) = (l-m 2 )(m 2 -2). (65) 

For E — > 0, we have (1 — m 2 ) sinh(a;o/3^/2) — * 8, a(m)E(m) + b(m)K(m) — > 2 and go;9i ~~ ► 
0, g 2 — > 1, so that A(/3) — > A\Ju)l/ixh, which reproduces the zero-temperature decay rate T . Let us 
turn now to the limit E — > 1, where the period behaves as 



2tt 



(3h= — 1 + — (1 — E) + 



36 



(66) 



The leading term gives a crossover temperature = huo/2n |17| , i.e. for > the decay 
rate is given by the familiar Arrhenuis-Kramers formula fi~8| |. On the other hand, for < the 
macroscopic tunneling through the barrier becomes more probable, and the decay rate is given by 
(|62|). Recalling the energy unit U defined by (f£9[), and (|2|) we find 



For small (f3 - f3 c ) / (3 C > 0, from 



2tt 
Tj and 



Uo 
k B 



5 1/4 (1 + 0(5 1/2 )). 



S, 



2nkB 

we obtain the bounce action 

2 



cl 



h 



P~^Pc 



P~Pc 
, Pc 



(67) 



(6£ 
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and 



77(P-p c \ 20867 



15hn 



\ 2 / 



20 



+ 



2400 



Pc 



(69) 




0.001 0.01 



100 1000 



Figure 6: The decay rate To as a function of the asymmetry parameter A for v = 953 sec , ao/a 
-2.13 x 10 3 and 5 = 5.0 x 10 -3 . 



4 Conclusion 

In this paper we have investigated the macroscopic tunneling of the metastable condensate of 7 Li. 
When the number of particles in the condensate exceeds a critical value N* = 



7i/2P*ao/a, the 

metastable condensate no longer exists and the equation giving critical points has no solutions. In a 
region extremely close to N*, i.e. 5 = 1 — P/P* <C 1, we have shown that the action takes a rather 



simple form (|5lD, and explicitly calculated the decay rate of the metastable condensate using the 
WKB approximation. 

Finally we make some remarks on our results. In order to justify the WKB approximation, we 
should choose the effective Plank constant h to satisfy the condition h -C 1. On the other hand, 
for very small h, it is impossible to observe the macroscopic tunneling; the formula ( |54] ) provides an 
estimate of the tunneling decay rate, To ~ 0{e~ l ^ h ). This implies rather severe conditions on the 
parameter 5 through the equation (0). If we use the experimental data at Rice University for the 

-2.13 x 10 3 and v 



trapping potential [|Tl], [12 
are given by h 



A ~ 0.867, a /a ~ —2.13 x 10 3 and v ~ 953 sec x , then the conditions 
2.92xilr 4 cT 5/4 < 1 and T ~ 8.17 x 10 5 5 7/8 exp (-6.72 x 10 3 5 5/4 ) ~ 0(1) sec" 1 . 
Consequently, we have a typical region 3.0 x 10~ 3 < 5 < 7.0 x 10~ 3 . Temperature effects on the 
tunneling decay rate are estimated by using the equations ( |68|) and (|69|): T(/3) is monotone decreasing 
for (3 > p c and hence Ar = r(/3) — r < r(/3 c ) — T . For instance, for S = 5.0 x 10 -3 the decay 
rate at zero temperature is r — 1.03 sec -1 and Ar < 2.79 sec -1 . The crossover temperature is then 
given by P~ l ~ 1.02 nK, which may be a realizable temperature in the experiments. The details of 
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a crossover region have been discussed in IDJ; there is a narrow crossover region of 0(h 3/2 ), where 
the decay rate is given by 



r(/3)T 



' 8u 3 . (Lu {3h\ 
— - — - sinn erf 

15/iTT 2 V 2 / 



/36 

Wc 



(P~Pc 



cxp 



-p+ 



5 I P e . 



(70) 



with erf(x) = (2rr) / dy exp(—y 2 /2). For very small h <C 10 2 , this formula matches smoothly 

J — oo 

onto (|62"D and T(0)T o = (Lj /2n) [smh(u ph/2)/ sm(uj ph/2)} exp(— (3) (Arrhenuis - Kramers formula) 
near P c . However, we can not apply the formula to the macroscopic tunneling since the value of h 
in our situation is too large. We leave the issue of crossover region for future research. The shape 
of the trapping potential also has some effect on the behavior of the decay rate IV as shown in 
Fig.|| the effect is significant for the disk-shaped potential (A > 1), although it is rather small for 
the cigar-shaped potential (A <C 1) and r is of order 10~ 3 sec -1 , independently of A. 
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